#ifndef ESDIRK_HPP
#define ESDIRK_HPP

#include <iostream>
#include <vector>
#include <math.h>
#include <fstream>
#include "norm.hpp"
#include "Function.hpp"

using namespace std;

class ESDIRK : public TimeIntegrator
{
private:
    int n_steps;
    double T, mu;
    vector<double> value;
    vector<vector<double>> result;
    Norm norm;

public:
    ESDIRK() {}
    ESDIRK(int n_, double time, vector<double> Val) : n_steps(n_), T(time), value(Val) {}
    ~ESDIRK() {}
    void solver();
    void error();
};

void ESDIRK::solver()
{
    auto epsilon = 0.000000001;
    ofstream outfile;
    outfile.open("ouput.m");
    outfile << "points = [";
    double k = 1.0 * T / n_steps;
    vector<double> U = value;
    outfile << U << endl;
    result.push_back(U);

    // result.push_back(U);
    for (int i = 0; i < n_steps; i++)
    {
        auto y1 = f(U, i * k);
        auto y2 = y1;
        while (true)
        {
            auto temp = f(U + (0.25 * k) * y1 + (0.25 * k) * y2, i * k + 0.5 * k);
            if (norm.inftynorm(y2 - temp) < epsilon)
                break;
            else
                y2 = temp;
        }
        auto y3 = y2;
        while (true)
        {
            auto temp = f(U + (8611 * k / 62500) * y1 - (1743 * k / 31250) * y2 + (0.25 * k) * y3, i * k + (83 * k / 250));
            if (norm.inftynorm(y3 - temp) < epsilon)
                break;
            else
                y3 = temp;
        }
        auto y4 = y3;
        while (true)
        {
            auto temp = f(U + (5012029 * k / 34652500) * y1 - (654441 * k / 2922500) * y2 + (174375 * k / 388108) * y3 + (0.25 * k) * y4, i * k + (31 * k / 50));
            if (norm.inftynorm(y4 - temp) < epsilon)
                break;
            else
                y4 = temp;
        }
        auto y5 = y4;
        while (true)
        {
            auto temp = f(U + (15267082809 * k / 155376265600) * y1 - (71443401 * k / 120774400) * y2 + (730878875 * k / 902184768) * y3 + (2285395 * k / 8070912) * y4 + (0.25 * k) * y5, i * k + (17 * k / 20));
            if (norm.inftynorm(y5 - temp) < epsilon)
                break;
            else
                y5 = temp;
        }
        auto y6 = y5;
        while (true)
        {
            auto temp = f(U + (82889 * k / 524892) * y1 + (15625 * k / 83664) * y3 + (69875 * k / 102672) * y4 - (2260 * k / 8211) * y5 + (0.25 * k) * y6, i * k + k);
            if (norm.inftynorm(y6 - temp) < epsilon)
                break;
            else
                y6 = temp;
        }
        U = U + (82889 * k / 524892) * y1 + (15625 * k / 83664) * y3 + (69875 * k / 102672) * y4 - (2260 * k / 8211) * y5 + (0.25 * k) * y6;
        outfile << U << endl;
        result.push_back(U);
    }
    outfile << "];" << endl;
    outfile << "scatter(points(:,1), points(:,2),\".\");" << endl;
    outfile.close();
}

void ESDIRK::error()
{
    // cout << result[n_steps] << endl;
    // cout << result[0] << endl;
    double error = norm.Inftynorm(result[n_steps] - result[0]);
    cout << error << endl;
}

TimeIntegrator* createESDIRK(int n_, int p_, double T_, vector<double> value_) { return new ESDIRK(n_, T_, value_); }

#endif